
* =============================================================================
* File Name: 2016_EC_EnvRegulation--RobustnessChecks_March2019.do

* File Description: This file contains the code for all robustness tests. 

* Need each robustness test for Z, X, exports, and net-exit. May also want to show results 
* for e, and potentially by productivity level?

* The identifiers are:
* 	- Plant identifier is NPRI_ID
* 	- Pollutant identifier is CAS_Number
* 	- Year identifier is REP_PERIOD 

* Creation date: March 6, 2019

* This version: March 6, 2019

* Author: Nouri Najjar
* =============================================================================

* Specify Pollutant
local x = "PM25"

* -----------------------------------------------------------------------------
* Section 1: Call the npri-asm data with air quality from the cder server. 

* This calls the npri-asm data (with CMA-year level air quality) that has been
* prepared for the analysis. To save memory space, all non-polluters are dropped.
* -----------------------------------------------------------------------------
* Set server folder. 
local datadir \\f4cder01\2016_EC_EnvRegulation\DATA\

* Open dataset. Use the ASM CMA Matched Data. 
use "`datadir'asm_npri_2002to2012`x'_withAQ_analysis.dta", clear
* -----------------------------------------------------------------------------

* -----------------------------------------------------------------------------
* Section 2: Prepare data for Dynamic DDD robustness test
* -----------------------------------------------------------------------------
*******
* Restrict Sample to Pollutant of Interest
*******
gen PM25=(cas_number=="NA - M10")
gen NOx=(cas_number=="11104-93-1")

keep if `x'==1
*******	 

***
* Create variable that tracks the number of years before/after the plant is treated. 
***
* First year plant is treated
bysort NPRI_ID: egen FirstYear_PM25Treated = min(cond(PM25Treated==1,REP_PERIOD,.))
* Years pre/post treatment 
gen YearsSince_PM25Treated = REP_PERIOD - FirstYear_PM25Treated
***

***
* Create indicator for number of years before treatment. 
***
forvalues i = 1/9 {
	gen YearsSince_PM25_Pre`i' = cond(YearsSince_PM25Treated==-`i',1,0)
}
***

***
* Create variables that track the number of years a plant is treated by 
* each standard. 
***
* Count number of years treated - PM
bysort NPRI_ID (REP_PERIOD): gen YearsTreated_PM25store = sum(PM25Treated) if PM25Treated==1 & dupPlantYear<=1
bysort NPRI_ID REP_PERIOD (YearsTreated_PM25store): gen YearsTreated_PM25=YearsTreated_PM25store[1]

* Create indicator that tracks the number of years in treatment. 
forvalues i = 1/10 {
	gen YearsTreated_PM25_`i' = cond(YearsTreated_PM25==`i',1,0)
}

* Group years after 4 years treated together
gen YearsTreated_PM25_GT4 = cond(YearsTreated_PM25>=4 & missing(YearsTreated_PM25)==0,1,0)
***


***
* Create variable "zero" to store the number zero for graphs. This is used in 
* both the dynamic DDD and flexible DDD robustness tests.
***
gen zero = 0 
***
* -----------------------------------------------------------------------------



* -----------------------------------------------------------------------------
* Section 3: Trim sample and create new variables  
* -----------------------------------------------------------------------------

*******
* Restrict Sample to 2004-2010
*******
keep if yr4>=2004 & yr4<=2010
*******

*******
* Create new industry variable that takes NAICS3 for all industries except
* 4-digit treated industries.
gen Industry = naics3
replace Industry = naics4 if naics4>= 3270 & naics4<=3279
replace Industry = naics4 if naics4>= 3240 & naics4<=3249

* Generate new industry year indicator
egen Ind3Year = group(Industry yr4)
*******

*******
* Drop all plants that either change cma over the sample or change industry
*******
bysort NPRI_ID: egen Mean_CMA=mean(CMA_ID_Merge)
gen ChangeCMA = cond(Mean_CMA!=CMA_ID_Merge,1,0)

bysort NPRI_ID: egen Mean_Industry=mean(Industry)
gen ChangeIndustry = cond(Mean_Industry!=Industry,1,0)

* Drop if change cma
drop if ChangeCMA==1
drop ChangeCMA
* drop if change industry
drop if ChangeIndustry==1
drop ChangeIndustry
*******


*******
* Create common sample across all regressions
*******
* Drop non-zero pollution
drop if missing(LN_total_air_converted)==1

* Drop all plants with missing sales or employment
drop if missing(LN_real_vsm)==1
drop if missing(LN_totalemp)==1

* Drop plants with missing VA or VA<0
drop if missing(LN_real_vam)==1 
drop if real_vam<0

* Drop all plants with missing CMA information
drop if missing(CMA_ID_Merge)==1

* Drop any plants with multiple ASM observations in a given year. 
bysort NPRI_ID cas_number yr4: gen dupLocNo = cond(_N==1,0,_n)
drop if dupLocNo>1
*******

*******
* Drop CMAs with 1 industry and industries with 1 cma
*******
* Tag distinct values
egen tag = tag(CMA_ID_Merge Industry)

* Drop if only one industry in CMA
egen distinct_CMA = total(tag), by(CMA_ID_Merge)
drop if distinct_CMA==1

* Drop if only one CMA in industry
egen distinct_Industry = total(tag), by(Industry)
drop if distinct_Industry==1
drop tag 
*******

*******
* Drop plants with only one year for each pollutant
*******
bysort NPRI_ID cas_number: egen NumberYears_PolPlant = count(NPRI_ID)
drop if NumberYears_PolPlant==1
*******

*******	
* CMA-Industry indicator for clustering standard errors
*******	
egen CMA_Ind3 = group(CMA_ID_Merge Industry)
*******	

*******	 
* Recode missing treatment plants as control
*******	 
recode PM25Treated (.=0)
recode O3Treated (.=0)
*******	 

* Change directory to save graphs. 
local datadir \\f4cder01\2016_EC_EnvRegulation\

* Keep regions with AQ data
keep if missing(PM25_CMA)==0
* -----------------------------------------------------------------------------



* -------------------------------------------------------------------------
* 1. Dynamic DDD robustness test
* -------------------------------------------------------------------------

***
* Regressions - PM emitters only. 
* Note that this also produces figures of the regression results. 
***
*Drop T-1
reghdfe LN_total_air_converted YearsSince_PM25_Pre3 YearsSince_PM25_Pre2 YearsTreated_PM25_1 YearsTreated_PM25_2 YearsTreated_PM25_3 YearsTreated_PM25_GT4 zero O3Treated [pw=Weights_`x'], absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'Pollution_dynamicDDD
reghdfe LN_real_vsm YearsSince_PM25_Pre3 YearsSince_PM25_Pre2 YearsTreated_PM25_1 YearsTreated_PM25_2 YearsTreated_PM25_3 YearsTreated_PM25_GT4 zero O3Treated [pw=Weights_`x'], absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'Sales_dynamicDDD
reghdfe LN_real_vam YearsSince_PM25_Pre3 YearsSince_PM25_Pre2 YearsTreated_PM25_1 YearsTreated_PM25_2 YearsTreated_PM25_3 YearsTreated_PM25_GT4 zero O3Treated [pw=Weights_`x'], absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'VA_dynamicDDD

* Print results in table
esttab `x'Pollution_dynamicDDD `x'Sales_dynamicDDD `x'VA_dynamicDDD using 2016_EC_EnvRegulation--`x'DynamicDDD_Table.xlsx, csv ///
	   keep(YearsSince_PM25_Pre3 YearsSince_PM25_Pre2 YearsTreated_PM25_1 YearsTreated_PM25_2 YearsTreated_PM25_3 YearsTreated_PM25_GT4) b(%9.3f) se(%9.3f)  r2 obslast star(* 0.10 ** 0.05 *** 0.01) coeflabels(YearsSince_PM25_Pre3 "{&le}T-3"  YearsSince_PM25_Pre2 "T-2" YearsTreated_PM25_1 "T" YearsTreated_PM25_2 "T+1" YearsTreated_PM25_3 "T+2" YearsTreated_PM25_GT4 "{&ge}  T+3") ///
	   mtitles("Emissions" "Sales" "Value Added") nonotes addnotes("Notes: Standard errors clustered by CMA-industry.") fragment replace 	  


* Graph results
local datadir \\f4cder01\2016_EC_EnvRegulation\

coefplot (`x'Pollution_dynamicDDD, label(Pollution)) (`x'Sales_dynamicDDD, label(Sales)) (`x'VA_dynamicDDD, label(Value Added)), level(90) omitted vertical yline(0, lcolor(red)) graphregion(color(white)) ytitle("Percent Change") xtitle("Years Pre/Post Treatment") xline(3.5, lcolor(red) lwidth(medthick) lpattern(dash)) ///
	keep(YearsSince_PM25_Pre3 YearsSince_PM25_Pre2 zero YearsTreated_PM25_1 YearsTreated_PM25_2 YearsTreated_PM25_3 YearsTreated_PM25_GT4) ///
	order(YearsSince_PM25_Pre3 YearsSince_PM25_Pre2 zero YearsTreated_PM25_1 YearsTreated_PM25_2 YearsTreated_PM25_3 YearsTreated_PM25_GT4) ///
	coeflabels(YearsSince_PM25_Pre3 = "{&le}T-3"  YearsSince_PM25_Pre2 = "T-2" zero = "T-1" YearsTreated_PM25_1 = "T" YearsTreated_PM25_2 = "T+1" YearsTreated_PM25_3="T+2" YearsTreated_PM25_GT4 = "{&ge}  T+3")
graph save "`datadir'Batchsubmit\2016_EC_EnvRegulation--`x'DynamicDDD_Graph.gph", replace
***

* Save data
keep if e(sample)
local datadir \\f4cder01\2016_EC_EnvRegulation\DATA\
saveold "`datadir'2016_EC_EnvRegulation--`x'MainRobustnessALL_March2019.dta", replace

* -------------------------------------------------------------------------


* -------------------------------------------------------------------------
* 2. Flexible triple difference robustness test

* Allows treatment to vary by city-year's AQ. 
* -------------------------------------------------------------------------
*****
* Create variables for regressions
*****

***
* Bins of PM air quality
***
gen PM25_15 = cond(PM25_CMA>=15 & PM25_CMA<18 | missing(PM25_CMA)==1,1,0)
gen PM25_18 = cond(PM25_CMA>=18 & PM25_CMA <21 & missing(PM25_CMA)==0,1,0)
gen PM25_21 = cond(PM25_CMA>=21 & PM25_CMA <24 & missing(PM25_CMA)==0,1,0)
gen PM25_24 = cond(PM25_CMA>=24 & PM25_CMA <27 & missing(PM25_CMA)==0,1,0)
gen PM25_26 = cond(PM25_CMA>=26 & PM25_CMA <28 & missing(PM25_CMA)==0,1,0)
gen PM25_28 = cond(PM25_CMA>=28 & PM25_CMA <30 & missing(PM25_CMA)==0,1,0)
gen PM25_30 = cond(PM25_CMA>=30 & PM25_CMA <32 & missing(PM25_CMA)==0,1,0)
gen PM25_32 = cond(PM25_CMA>=32 & PM25_CMA <34 & missing(PM25_CMA)==0,1,0)
gen PM25_34 = cond(PM25_CMA>=34 & missing(PM25_CMA)==0,1,0)

***
***
* Bins of AQ interacted with target ind
***
gen PM25Treated_bin15 = PM25_15*TargetInd
gen PM25Treated_bin18 = PM25_18*TargetInd
gen PM25Treated_bin21 = PM25_21*TargetInd
gen PM25Treated_bin24 = PM25_24*TargetInd
gen PM25Treated_bin26 = PM25_26*TargetInd
gen PM25Treated_bin28 = PM25_28*TargetInd
gen PM25Treated_bin30 = PM25_30*TargetInd
gen PM25Treated_bin32 = PM25_32*TargetInd
gen PM25Treated_bin34 = PM25_34*TargetInd
***

***
* PM Regressions
* Also produces figures of regression results. 
***
reghdfe LN_total_air_converted PM25Treated_bin24 PM25Treated_bin26 PM25Treated_bin28 PM25Treated_bin30 PM25Treated_bin32 PM25Treated_bin34 zero [pw=Weights_`x'] if PM25_CMA>=20 & PM25_CMA<=50, absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'Pollution_flexibleDDD
reghdfe LN_real_vsm PM25Treated_bin24 PM25Treated_bin26 PM25Treated_bin28 PM25Treated_bin30 PM25Treated_bin32 PM25Treated_bin34 zero [pw=Weights_`x'] if PM25_CMA>=20 & PM25_CMA<=50, absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'Sales_flexibleDDD
reghdfe LN_real_vam PM25Treated_bin24 PM25Treated_bin26 PM25Treated_bin28 PM25Treated_bin30 PM25Treated_bin32 PM25Treated_bin34 zero [pw=Weights_`x'] if PM25_CMA>=20 & PM25_CMA<=50, absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'VA_flexibleDDD

* Print results in table
esttab `x'Pollution_flexibleDDD `x'Sales_flexibleDDD `x'VA_flexibleDDD using 2016_EC_EnvRegulation--`x'FlexibleDDD_Table.xlsx, csv ///
	   keep(PM25Treated_bin24 PM25Treated_bin26 PM25Treated_bin28 PM25Treated_bin30 PM25Treated_bin32 PM25Treated_bin34) b(%9.3f) se(%9.3f)  r2 obslast star(* 0.10 ** 0.05 *** 0.01) coeflabels(PM25Treated_bin24 "24-25" PM25Treated_bin26 "26-27" PM25Treated_bin28 "28-29" PM25Treated_bin30 "30-31" PM25Treated_bin32 "32-33" PM25Treated_bin34 "{&ge}34") ///
	   mtitles("Emissions" "Sales" "Value Added") nonotes addnotes("Notes: Standard errors clustered by CMA-industry.") fragment replace 	  

* Print results as graph
local datadir \\f4cder01\2016_EC_EnvRegulation\

coefplot (`x'Pollution_flexibleDDD, label(Pollution)) (`x'Sales_flexibleDDD, label(Sales)) (`x'VA_flexibleDDD, label(Value Added)), level(90) omitted vertical yscale(range(-0.1 0.5)) yline(0, lcolor(red)) graphregion(color(white)) ytitle("Percent Change") xtitle("CMA Air Quality") ///
	keep(zero PM25Treated_bin24 PM25Treated_bin26 PM25Treated_bin28 PM25Treated_bin30 PM25Treated_bin32 PM25Treated_bin34) ///
	order(zero PM25Treated_bin24 PM25Treated_bin26 PM25Treated_bin28 PM25Treated_bin30 PM25Treated_bin32 PM25Treated_bin34) ///
	coeflabels(zero = "{&le}20" PM25Treated_bin24 = "24-25" PM25Treated_bin26 = "26-27" PM25Treated_bin28 = "28-29" PM25Treated_bin30 = "30-31" PM25Treated_bin32 = "32-33" PM25Treated_bin34 = "{&ge}34") ///
	addplot(scatteri -1 3.5  -1 4.5, recast(area) fcolor(orange*0.2) lcolor(orange*0.2) below || scatteri 0.5 3.5  0.5 4.55, recast(area) fcolor(orange*0.2) lcolor(orange*0.2) below || scatteri -1 4.5  -1 7.5, recast(area) fcolor(red*0.4) lcolor(red*0.4) below || scatteri 0.5 4.5  0.5 7.5, recast(area) fcolor(red*0.4) lcolor(red*0.4) below) ///
	text(0.3 4 "Warning" "Range" 0.3 6.1 "Violation Range")
graph save "`datadir'Batchsubmit\2016_EC_EnvRegulation--`x'FlexibleDDD_Graph.gph", replace
***
* -------------------------------------------------------



/* OMITTED FROM RESULTS
* -------------------------------------------------------------------------
* 3. Estimation as separate Diff in Diffs for targeted and non-targeted industries. 
* Standard errors are clustered by CMA

* Variables
*	PM25_AboveThreshold = indicator for CMA-years violating PM standard (defined in other file)
* 	TargetIndv2 = indicator for targeted industries (defined in other file)
* 	CMAUID = CMA ID
* -------------------------------------------------------------------------

* Target Ind only
reghdfe LN_total_air_converted PM25_AboveThreshold [pw=Weights_`x'] if TargetInd==1, absorb(NPRI_ID yr4) vce(cluster CMA_ID_Merge)
estimates store `x'Emissions_DD_Target
reghdfe LN_real_vsm PM25_AboveThreshold [pw=Weights_`x'] if TargetInd==1, absorb(NPRI_ID yr4) vce(cluster CMA_ID_Merge)
estimates store `x'Sales_DD_Target
reghdfe LN_real_vam PM25_AboveThreshold [pw=Weights_`x'] if TargetInd==1, absorb(NPRI_ID yr4) vce(cluster CMA_ID_Merge)
estimates store `x'VA_DD_Target
reghdfe LN_totalemp PM25_AboveThreshold [pw=Weights_`x'] if TargetInd==1, absorb(NPRI_ID yr4) vce(cluster CMA_ID_Merge)
estimates store `x'Emp_DD_Target

* Non Target Industries
reghdfe LN_total_air_converted PM25_AboveThreshold [pw=Weights_`x'] if TargetInd==0, absorb(NPRI_ID yr4) vce(cluster CMA_ID_Merge)
estimates store `x'Emissions_DD_NonTarget
reghdfe LN_real_vsm PM25_AboveThreshold [pw=Weights_`x'] if TargetInd==0, absorb(NPRI_ID yr4) vce(cluster CMA_ID_Merge)
estimates store `x'Sales_DD_NonTarget
reghdfe LN_real_vam PM25_AboveThreshold [pw=Weights_`x'] if TargetInd==0, absorb(NPRI_ID yr4) vce(cluster CMA_ID_Merge)
estimates store `x'VA_DD_NonTarget
reghdfe LN_totalemp PM25_AboveThreshold [pw=Weights_`x'] if TargetInd==0, absorb(NPRI_ID yr4) vce(cluster CMA_ID_Merge)
estimates store `x'Emp_DD_NonTarget


esttab `x'Emissions_DD_Target `x'Sales_DD_Target `x'VA_DD_Target `x'Emissions_DD_NonTarget `x'Sales_DD_NonTarget `x'VA_DD_NonTarget using 2016_EC_EnvRegulation--`x'SplitDDD_March2019.xlsx, tex ///
	   keep(PM25_AboveThreshold) b(%9.3f) se(%9.3f)  r2 obslast star(* 0.10 ** 0.05 *** 0.01) coeflabel(PM25_AboveThreshold "Violate PM25") ///
	   mtitles("Emissions" "Sales" "VA" "Emissions" "Sales" "VA") nonotes addnotes("Notes: Standard errors clustered by CMA-industry.") fragment replace 	  
* -------------------------------------------------------------------------
*/




   
* -------------------------------------------------------------------------
* 5. Allow for separate CMA-Year-FEs for different size emitters. 
* -------------------------------------------------------------------------
****
* Prepare dataset
****
* Total releases by CMA Year
bysort CMA_ID_Merge yr4: egen TotPlantPM25Rel_City = total(total_air_converted)

* Fraction of total city releases each plant accounts for
gen FracPM25 = total_air_converted/TotPlantPM25Rel_City

* Bins for fixed effect 
gen FracPM25Bin = 1 if cond(FracPM25<.01,1,0)
replace FracPM25Bin = 2 if cond(FracPM25>.01 & FracPM25<.20,1,0)
replace FracPM25Bin = 3 if cond(FracPM25>.20 & missing(FracPM25)==0,1,0)

* Fixed effect
egen FracPM25Bin_CMA_Year = group(FracPM25Bin CMA_Year)
****


****
* Estimates for paper
****
preserve
* Pollution
reghdfe LN_total_air_converted PM25Treated [pw=Weights_`x'], absorb(NPRI_ID Ind3Year FracPM25Bin_CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'Pol_BigEmitCMAYearFEs
reghdfe LN_real_vsm PM25Treated [pw=Weights_`x'], absorb(NPRI_ID Ind3Year FracPM25Bin_CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'Sales_BigEmitCMAYearFEs
reghdfe LN_real_vam PM25Treated [pw=Weights_`x'], absorb(NPRI_ID Ind3Year FracPM25Bin_CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'VA_BigEmitCMAYearFEs

* Save data
keep if e(sample)
local datadir \\f4cder01\2016_EC_EnvRegulation\DATA\
saveold "`datadir'2016_EC_EnvRegulation--`x'MainRobustnessTrend_March2019.dta", replace
restore

* Print all results in a tex file. 
esttab `x'Pol_BigEmitCMAYearFEs `x'Sales_BigEmitCMAYearFEs `x'VA_BigEmitCMAYearFEs using 2016_EC_EnvRegulation--`x'BigEmitterTrends_March2019.xlsx, csv ///
	   title("ln(Emissions)") keep(PM25Treated) b(%9.3f) se(%9.3f)  r2 obslast star(* 0.10 ** 0.05 *** 0.01) coeflabel(PM25Treated "PM 2.5 Standard") ///
	   mtitles("Emissions" "Sales" "Value Added") nonotes addnotes("Notes: Standard errors clustered by CMA-industry.") fragment replace 	  
* -------------------------------------------------------------------------



* -------------------------------------------------------------------------
* 6. Drop large emitters. Uses variables created in robustness check 5.
* -------------------------------------------------------------------------
* Estimates - drop above 20%
preserve
reghdfe LN_total_air_converted PM25Treated [pw=Weights_`x'] if FracPM25<=.2, absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'Emissions_Below20
reghdfe LN_real_vsm PM25Treated [pw=Weights_`x'] if FracPM25<=.2, absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'Sales_Below20
reghdfe LN_real_vam PM25Treated [pw=Weights_`x'] if FracPM25<=.2, absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'VA_Below20
keep if e(sample)
local datadir \\f4cder01\2016_EC_EnvRegulation\DATA\
saveold "`datadir'2016_EC_EnvRegulation--`x'DropEmitters_Col1_March2019.dta", replace
restore
* Estimates - drop above 10%
preserve
reghdfe LN_total_air_converted PM25Treated [pw=Weights_`x'] if FracPM25<=.1, absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'Emissions_Below10
reghdfe LN_real_vsm PM25Treated [pw=Weights_`x'] if FracPM25<=.1, absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'Sales_Below10
reghdfe LN_real_vam PM25Treated [pw=Weights_`x'] if FracPM25<=.1, absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'VA_Below10
keep if e(sample)
local datadir \\f4cder01\2016_EC_EnvRegulation\DATA\
saveold "`datadir'2016_EC_EnvRegulation--`x'DropEmitters_Col2_March2019.dta", replace
restore
* Estimates - drop above 5%
preserve
reghdfe LN_total_air_converted PM25Treated [pw=Weights_`x'] if FracPM25<=.05, absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'Emissions_Below5
reghdfe LN_real_vsm PM25Treated [pw=Weights_`x'] if FracPM25<=.05, absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'Sales_Below5
reghdfe LN_real_vam PM25Treated [pw=Weights_`x'] if FracPM25<=.05, absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'VA_Below5
keep if e(sample)
local datadir \\f4cder01\2016_EC_EnvRegulation\DATA\
saveold "`datadir'2016_EC_EnvRegulation--`x'DropEmitters_Col3_March2019.dta", replace
restore
* Estimates - drop above 1%
preserve
reghdfe LN_total_air_converted PM25Treated [pw=Weights_`x'] if FracPM25<=.01, absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'Emissions_Below1
reghdfe LN_real_vsm PM25Treated [pw=Weights_`x'] if FracPM25<=.01, absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'Sales_Below1
reghdfe LN_real_vam PM25Treated [pw=Weights_`x'] if FracPM25<=.01, absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'VA_Below1
keep if e(sample)
local datadir \\f4cder01\2016_EC_EnvRegulation\DATA\
saveold "`datadir'2016_EC_EnvRegulation--`x'DropEmitters_Col4_March2019.dta", replace
restore

* Print all results in a tex file. 
esttab `x'Emissions_Below20 `x'Sales_Below20 `x'VA_Below20 `x'Emissions_Below10 `x'Sales_Below10 `x'VA_Below10 `x'Emissions_Below5 `x'Sales_Below5 `x'VA_Below5 `x'Emissions_Below1 `x'Sales_Below1 `x'VA_Below1 using 2016_EC_EnvRegulation--`x'DropLargeEmitters_March2019.xlsx, csv ///
	   title("ln(Emissions)") keep(PM25Treated) b(%9.3f) se(%9.3f)  r2 obslast star(* 0.10 ** 0.05 *** 0.01) coeflabel(PM25Treated "PM 2.5 Standard") ///
	   mtitles("Emissions" "Sales" "Value Added" "Emissions" "Sales" "Value Added" "Emissions" "Sales" "Value Added" "Emissions" "Sales" "Value Added") nonotes addnotes("Notes: Standard errors clustered by CMA-industry.") fragment replace 	  
* ---------------------------------------------------------------------


* -------------------------------------------------------------------------
* 7. Allow treatment effect to vary for multiplant and single plant firms
* -------------------------------------------------------------------------
****
* Create variables for multiplant tests
****
* Count number of plants in same enterprise
bysort meadent yr4: egen NumberPlants = count(yr4)

* Indicator for mulitplant firms
gen Multiplant = cond(NumberPlants>1,1,0)
tab Multiplant

* Interact treatment with multiplant
gen PM25Treated_Multi = PM25Treated*Multiplant
****

****
* Estimates for paper
****
preserve
* Emissions
reghdfe LN_total_air_converted PM25Treated_Multi PM25Treated Multiplant [pw=Weights_`x'], absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'Emissions_Multi
* Sales
reghdfe LN_real_vsm PM25Treated_Multi PM25Treated Multiplant [pw=Weights_`x'], absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'Sales_Multi
* Value Added
reghdfe LN_real_vam PM25Treated_Multi PM25Treated Multiplant [pw=Weights_`x'], absorb(NPRI_ID Ind3Year CMA_Year) vce(cluster CMA_Ind3)
estimates store `x'VA_Multi

* Save data
keep if e(sample)
local datadir \\f4cder01\2016_EC_EnvRegulation\DATA\
saveold "`datadir'2016_EC_EnvRegulation--`x'MainRobustness_Multiplant_March2019.dta", replace
restore

* Print all results in a tex file. 
esttab `x'Emissions_Multi `x'Sales_Multi `x'VA_Multi using 2016_EC_EnvRegulation--`x'MultiplantChecks_March2019.xlsx, csv ///
	   title("ln(Emissions)") keep(PM25Treated PM25Treated_Multi Multiplant) order(PM25Treated PM25Treated_Multi Multiplant) b(%9.3f) se(%9.3f)  r2 obslast star(* 0.10 ** 0.05 *** 0.01) coeflabel(PM25Treated "PM 2.5 Standard" PM25Treated_Multi "PM x Multi" Multiplant "Multiplant") ///
	   mtitles("Emissions" "Sales" "Value Added") nonotes addnotes("Notes: Standard errors clustered by CMA-industry.") fragment replace 	  
* -------------------------------------------------------------------------


* -----------------------------------------------------------------------------



